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Abstract 

We present four-dimensional ab initio potential energy surfaces for the three different spin states 
of the NH( 3 £~) - NH( 3 S~) complex. The potentials are partially based on the work of Dhont 
et al. [J. Chem. Phys. 123, 184302 (2005)]. The surface for the quintet state is obtained at the 
RCCSD(T)/aug-cc-pVTZ level of theory and the energy differences with the singlet and triplet 
states are calculated at the CASPTn/aug-cc-pVTZ (n = 2,3) level of theory. The ab initio 
potentials are fitted to coupled spherical harmonics in the angular coordinates, and the long range 
is further expanded as a power series in 1/R. The RCCSD(T) potential is corrected for a size- 
consistency error of about 0.5 x 10" 6 E h prior to fitting. The long-range coefficients obtained 
from the fit are found to be in good agreement with first and second-order perturbation theory 
calculations. 
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I. INTRODUCTION 



The field of cold (T < 1 K) and ultracold (< 1 mK) molecules has attracted great 
interest in the last few years. The production of such (ultra) cold species may find important 
applications in condensed matter physics jl|, high precision spectroscopy j^, 3, 4], physical 
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11] . There are, in principle, two 



chemistry 0, [(j, Q, g[ Q], and quantum computing 
different strategies for producing molecular samples at (ultra) low temperatures. In indirect 
methods, cold molecules are formed by pairing up atoms that are already cooled down to 
the ultracold regime. Examples of such methods include photoassociation Il2| and Feshbach 
association [lj|. Conversely, direct methods such as Stark deceleration and buffer gas 
cooling [ljj] employ a scheme in which pre-existing molecules are cooled down from higher 
temperatures. 

One of the most promising candidates for direct-cooling experiments is the NH radical. 
NH(X 3 £~) has a relatively large magnetic moment of 2 un, making it suitable for buffer 
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181 ] . Furthermore, the metastable 



gas cooling and magnetic trapping experiments 
a 1 A state of NH, which exhibits a linear Stark effect, can be efficiently Stark decelerated 
and trapped in an electrostatic field. Subsequent excitation of the A 3 U <— a 1 A transition 
followed by spontaneous emission to the ground state yields cold NH(X 3 £~) molecules, 
which in turn may be trapped in a mag netic field Q, Q. This scheme also allows for 
reloading of the magnetic trap, thus providing a means to increase phase-space density. 

At present, direct-cooling methods for NH are limited to temperatures of a few hundred 
mK. If the density of trapped molecules is sufficiently high, it may be possible to reach the 
ultracold regime by means of evaporative cooling. This process relies on elastic NH + NH 
collisions as the trap depth is gradually reduced. Inelastic spin-changing collisions between 
trapped NH molecules will lead to immediate trap loss and are therefore unfavorable. It is 
generally accepted that, in order for evaporative cooling to be successful, elastic collisions 
should be a few orders of magnitude more efficient than inelastic transitions 
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221 ] . In the case of NH(X 3 S~), the only magnetically trappable state is the low- field seeking 
Ms = 1 state, with Ms denoting the spin projection quantum number. A collision complex 
of two such molecules is in the Ms = 2 level of the NH-NH high-spin quintet (S = 2) state. 
Inelastic collisions between NH molecules may either change the Ms quantum number of the 
quintet state, or change the total spin S to produce singlet or triplet complexes. The S = 



2 



and 1 dimer states are chemically reactive 23j, |24| and, although unfavorable for evaporative 



cooling, could be of interest in the context of cold controlled chemistry [9|. 

A recent theoretical study by Kajita 25|, in which only the electric dipole-induced dipole 
and magnetic dipole-dipole interactions were considered, showed that evaporative cooling of 
NH is likely to be feasible. A more rigorous quantum calculation of elastic and inelastic cross 
sections, however, requires knowledge of the full NH-NH interaction potentials for all three 
spin states. In particular the long-range potential, which governs the dynamics at (ultra) low 
temperatures, should be described very accurately. For NH-NH the dominant long-range 
term is the electrostatic dipole-dipole interaction, which scales with the intermolecular dis- 
tance R as R~ 3 . If, however, the molecules are freely rotating, all multipole-multipole terms 
average out to zero and the isotropic (-R -6 ) dispersion and induction interactions become 
important. 

Dhont et al. have recently constructed four- dimensional ab initio potential energy 
surfaces for NH-NH which, in principle, contain all relevant long range contributions. They 
employed the partially spin-restricted coupled-cluster method with single and double excita- 
tions and a perturbative treatment of triples [RCCSD(T)] 27|, |28[ to obtain the surface for 
the NH-NH quintet state. We found, however, that this surface exhibits erroneous behavior 
in the long range due to a lack of size consistency in the open-shell RCCSD(T) method. In 
the present paper, we report more accurate ab initio calculations that are corrected for this 
undesirable feature, and which allow for an analytical fit of the long-range potential. The 
fit of the short-range potentials is also improved. 

This paper is organized as follows. In Section III Al we first address the RCCSD(T) 
size-consistency problem and present new RCCSD(T) calculations for the long range of the 
NH-NH potential. Long-range perturbation theory calculations are discussed in Section 
III B| and new CASPTn (n = 2, 3) calculations for the short range of the singlet and triplet 
potentials are presented in Section III CI The fit of the different potentials is described in 
Section UTTl followed by a discussion of the results in Section HVl Finally, conclusive remarks 
are given in Section |VJ 
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II. ELECTRONIC STRUCTURE CALCULATIONS 



A. RCCSD(T) potential energy surface 

The coupled-cluster (CC) approach is one of the most accurate ab initio methods available 
for calculating potential energy surfaces. This method requires a single Slater determinant 
function as the reference state, which in the case of NH-NH implies that only the high-spin 
quintet state is suitable for coupled-cluster calculations. At large intermolecular distances 
however, the energy splittings between the three different spin states become negligible, and 
thus the CC potential also applies to the singlet and triplet states at long range. In this 



section, we will show that the previously reported NH-NH RCCSD(T) potential [26J contains 
a size-consistency error that becomes apparent at large R. We also present new ab initio 
calculations that are corrected for this defect. The coordinates used to describe the NH-NH 
potential energy surfaces are the four intermolecular Jacobi coordinates (R, 6 a, Ob, <fi)- The 
coordinate R is the length of the intermolecular vector R that connects the centers of mass 
of monomers A and B, 6a and &b are the polar angles of the NH monomer axes relative to 
R, and (j) is the dihedral angle between the planes through R and the monomer axes (see 



also Fig. 1 of Ref. [26]). All interaction potentials are computed using the supermolecule 



approach with the counterpoise correction method of Boys and Bernardi 29]. 



1. Size consistency 

It is well established that coupled-cluster theory for closed-shell systems is rigorously 
size-consistent. For open-shell species, however, where the problem of nonzero spin arises, 
this issue is not straightforward. It was demonstrated in 2006 by Heckert et al. {30} that 
several spin-adapted CCSD schemes applied to the triplet F( 2 P) - F( 2 P) system exhibit 
size-consistency errors on the order of 10~ 7 - 10~ 8 Eh- The reason for this is still unclear, 
)ut it has been suggested that the problem lies in the truncation of the cluster operator 
3o| . Although the errors are very small, the effect becomes apparent when considering 
interactions at low temperatures, where the total energy of the system may be of a similar 
order of magnitude (10~ 7 Eh ~ 0.03 K). Hence, a lack of size consistency imposes a significant 
limitation on the accuracy of calculations in the (ultra)cold regime. 



When evaluating the NH( 3 £~) - NH( 3 £~) quintet potential of Ref. [26] in more detail, 



we indeed found that the interaction energy does not tend to zero at large intermolecular 
distances. At R — 30 000 a , the size-consistency error is —4.8823 x 10~ 6 Eh calculated 
at the RCCSD level of theory with the augmented correlation-consistent polarized valence 
triple-zeta (aug-cc-pVTZ) basis set |3l|, and +0.5129 x 10~ 6 E h at the RCCSD (T)/aug-cc- 
pVTZ level of theory. It should be noted that these errors are independent of the relative 
orientation of the monomers, i.e., the lack of size consistency affects only the isotropic part 
of the potential. The results for other basis sets are given in Table [H It can be seen that 
the error is largest at the RCCSD level and increases with the size of the basis set. The 
inclusion of triple excitations reduces the error by approximately one order of magnitude 
and, for most basis sets, also changes its sign. 

Although the problem has not been solved yet, we found that the NH-NH RCCSD (T) 
potential can be easily corrected for the lack of size consistency by simply subtracting the 
error, calculated at 30 000 do, from all ab initio points. We compared these corrected energies 
with the results obtained from a spin-unrestricted CCSD(T) [UCCSD(T)] calculation, of 
which the energies do converge to zero at long range [i.e. UCCSD(T) is size-consistent]. At 
R = 30.0 ao, the root-mean-square (RMS) difference between the UCCSD(T) and corrected 
RCCSD(T) data was calculated to be 9.1 xl0~ 9 Eh (0.08% of the mean absolute value of 
the potential) for a grid of 126 ab initio points. Without the size-consistency correction this 
difference would be 5.1 xl0~ 7 Eh (4.4%). Thus, the error subtraction at the RCCSD(T) 
level leads to significantly better agreement with the size-consistent UCCSD(T) method. 
Similar results were obtained at an intermolecular distance of 15.0 ao, where the RMS 
difference between the corrected RCCSD(T) and UCCSD(T) data is 7.0xl0~ 8 E h (0.07% of 
the mean absolute energy), as opposed to 5.4xl0~ 7 E h (0.54%) without the correction. 
At even smaller distances, the size-consistency error will become increasingly negligible 
compared to the total interaction energy, thus the correction will leave the short-range 
potential essentially unaffected. Based on these findings, we conclude that subtracting the 
error from all RCCSD (T) points does not significantly alter the accuracy of the potential, 
but does give the desired asymptotic behavior at long range. 
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2. Long-range RCCSD(T) calculations 



Although the size-consistency correction already constitutes an important refinement to 
the RCCSD(T) potential of Ref. 26(, we chose to improve the long range even further by 
performing new ab initio calculations. This is motivated by our aim to study collisions in the 
limit of zero temperature, for which it is desirable to have the long range in analytical form. 
In order to perform an accurate analytical fit, however, we found that the long-range ab initio 
energies should be converged to less than 10~ 10 Eh, while the data presented in Ref. [26( have 
been converged to only 10~ 8 E^. We therefore recalculated the points at large R with much 
tighter convergence thresholds, as low as 10 -13 Eh, to ensure that the fit will not be affected 
by numerical noise. The radial grid consisted of 8 points, approximately logarithmically 
spaced at 8.3, 10.0, 12.0, 14.4, 17.3, 20.8, 25.0, and 30.0 a . For the angular grid, we chose 
an 11-point Gauss- Legendre quadrature grid in (9a, 0b) an d an 11-point Gauss-Chebyshev 
grid in 0. These are known to be the most accurate quadratures on their respective domains 
321 ] . Due to the symmetry of the complex, only points with 9a + 9b < ^ arid < <fi < n were 



required in the calculations 



261 ] . The monomers were treated as rigid rotors, with the NH 



33|. The RCCSD(T) 



bond length fixed to the experimental equilibrium value of 1.0362 A 
energies were computed using the aug-cc-pVTZ basis set, with additional bond functions 
located at the midpoint of the intermolecular vector R (exponents s,p: 0.9, 0.3, and 0.1; 
d, /: 0.6 and 0.2; g: 0.3). All calculations were performed with the MOLPRO package 34]. 
As explained above, the size-consistency error of 0.51290xl0~ 6 Eh was subtracted from all 
RCCSD(T) points to ensure that the long range converges to zero. 



B. Perturbation theory calculations 

As an additional test for the accuracy of the RCCSD(T) long-range potential, we com- 
puted the long-range coefficients directly from first and second-order perturbation theory 



with the multipole expansion of the interaction operator (see e.g. Ref. [35j). The first-order 
(electrostatic) coefficients are expressed in terms of the permanent NH multipole moments, 
while the second-order (induction and dispersion) coefficients depend also on the static 
and dynamic polarizabilities of NH. The permanent multipole moments were obtained from 
finite field calculations at the RCCSD(T)/aug-cc-pVTZ level of theory and from density 
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functional theory (DFT), yielding two different sets of first-order coefficients. All DFT cal- 
culations were performed with the PBEO density functional j^f]] and the aug-cc-pVQZ basis 
set. The Kohn-Sham orbitals were obtained from a spin-restricted calculation using the 
DALTON program [37| • The Fermi- Amaldi asymptotic correction 38| was employed to im- 
prove the description of the NH densities. The ionization potentials used for this correction 
were taken from Ref. 3j|. For the static and dynamic NH polarizabilities, we performed 
spin-restricted time-dependent coupled Kohn-Sham (CKS) calculations 39]. Previous stud- 
ies have shown that CKS methods yield accurate van der Waals coefficients, comparable to 
the accuracies obtained with the best ab initio methods, for systems such as He2, Ne2, H 2 



dimer 



40] , and the open-shell O2 dimer 411 ] . The static polarizabilities and dynamic polar- 



izabilities at imaginary frequencies were obtained with a modified version of the SAPT2008 
package Q|, extended to treat open-shell fragments. Finally, the second-order long-range 
coefficients were computed from the DFT multipole moments and response functions using 
the POLCOR program [43] . 



C. CASPTn calculations 



As mentioned before, the RCCSD(T) quintet potential can also be used to describe the 
singlet and triplet NH-NH states at long range. In the short range, however, these lower- 
spin states must be treated with a different ab initio method. Dhont et al. 26j] employed the 
Complete Active Space with nth-order Perturbation Theory (CASPTn) method (n = 2, 3) 
to calculate the energy differences between the quintet state and the 5* = and 1 states, and 
added those to the RCCSD(T) quintet surface to obtain the singlet and triplet potentials: 



V s — V s — v s=2 

v n — v CASPTn v CASPTn 



, V S=2 
^ V RCCSD(T)' 



(1) 



When fitting the CASPTn energy splittings, which decay exponentially as a function of R, 
we found that the convergence thresholds used in Ref. 26j were not sufficiently stringent to 
reach the same accuracy as in the long range. Hence, we recalculated the CASPTn energies 
for all three spin states with much tighter convergence criteria. The active space consisted of 
the four orbitals that are singly occupied in the quintet state. The g 4 operator [44] was used 
to obtain size-consistent results, and a level shift of 0.4 was applied to enforce convergence. 
The interaction energies were computed for R = 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 
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9.0, 10.0, 12.0, and 14.4 a , with the energy threshold set to 10 13 E h for the points at 8.0 - 
14.4 a , 10" 12 E h at 7.0 and 7.5 a , 10~ u E h at 6.0 and 6.5 a , 10" 10 E h at 5.0 and 5.5 a , 
10 -9 Eh at 4.5 a , and 10 -8 Eh at 4.0 a . For the angular grid we used the same points as for 
the long-range RCCSD(T) calculations, i.e. an 11-point Gauss-Legendre quadrature in (9a, 
9 B ) and an 11-point Gauss- Chebyshev grid in 0. The CASPTn calculations were performed 
with MOLPRO [34] using the aug-cc-pVTZ basis set supplemented with bond functions. It 
should be noted that three points at 4.0 a failed to converge due to the strongly repulsive 
nature of the potential at small R. 



III. ANALYTICAL REPRESENTATION 

All three interaction potentials can be represented as follows: 

v{R,e A ,e B ,<f>) = v l a ^l(R)a La ,l b ,l(0a, 0b A) (2) 

L A ,L B ,L 

= V L a ,L b ,m( R ) A L a ,L b ,m(0a,0 B ,4>). (3) 

L A ,L B ,M 

The angular functions A LaLbL (9a, 9 b , 0) are defined as 

min( ^ Ls) I L A L B L\ 
Al a ,l b ,l(0a,0 b ,4>) = > v \C La ,m{Qa, 4>a)Cl b .~m(0 b , 4> b ), 

m=-±(l a ,l b ) \ M -M J 

min(L A ,L B ) / j j ^ \ 

E ^ A B \Al a ,l b A0a,0 b A), (4) 



A/=0 



M -M 



where C^^u(0, 0) are Racah-normalized spherical harmonics and = 0a — <Pb is the difference 
between the azimuthal angles of monomers A and B. The factor in brackets denotes a Wigner 
three- j symbol. The 'primitive' angular functions Al a ,l b ,m(0a,0b, (ft) are given by 

A La ,l b ,m(6a, b , 0) = Pl a ,m( cos a )Pl b ,m( cos 0b) cos M0, (5) 
where Pl,m (cos 9) are Schmidt semi- normalized associated Legendre functions defined in 



Ref. 



26) . The -R-dependent expansion coefficients are related to each other as 26] 



min(L A ,L B ) , ^ ^ j 

vl a ,l b ,l(R) = (2L + 1) (- 1 ) M ( 2 " M ( * * n ) "W^mW (6) 

M=0 



M — M 



S 



A. Long-range potential 



For the analytical long-range interaction, we use Eq. ([2]) and further expand the 
v l a ,l b ,l(R) coefficients in a power series in 1/R: 

v La , Lb ar) = Y, ~ CLA R t B ' L ' r ' " ■ (7) 

n 

Our choice of an 11-point Gauss-Legendre quadrature in (9 a, Ob) and an H-point Gauss- 
Chebyshev quadrature in ensures that the angular A La ^ Lb ^ l functions, when evaluated on 
the quadrature grid with the appropriate weights, are mutually orthogonal for all values of 
La and Lb up to 10 inclusive. Thus, we can perform the analytical fit in R [Eq. ([7])] for 
each (La, Lb, L) term separately. The values of n follow from a consideration of the possible 
first-order (electrostatic) and second-order (induction/dispersion) contributions (see e.g. Ref. 

for details). For the electrostatic terms, we have La + L B = L and n = La + L B + 1, 
with La > 1 and Lg > 1. The minimum value of 1 comes from the fact that the lowest 
nonvanishing permanent multipole moment of NH is the dipole. In the case of induction 
and dispersion interactions, La and Lb correspond to the order of two coupled multipole 
moments on monomers A and B, respectively. That is, La = \Ia — . . . ,1a + l' A and 
L b = \Ib — 1'b\, ■ ■ ■ + I'b, where I a, I 1 a, Ib, and l' B denote the orders of the uncoupled 
monomer multipole moments. La and Lb are in turn coupled to all possible L values, and 
for each (La,Lb,L) term we have n = Ia + I'a + + I'b + 2. Finally, due to the inversion 
symmetry of the total system, it can be shown that La + Lb + L is even, and since each 
monomer is a linear £ state molecule, I a + I'a + La and Ib + I'b + Lb must also be even 45]. 

The C La,Lb ,L,n fit coefficients of Eq. ([7j) were calculated as follows. For each set of 
(La, Lb, L) values, we first computed the lowest possible values of n in both first and second 
order. Since our long-range ab initio calculations were performed on a grid of eight R 
points, we could include a maximum of eight R~ n functions in the fit. We then fitted the 
size-consistency corrected RCCSD(T) data to the expansion of Eq. (T5]), and subsequently 
fitted each vl a ,l b ,l expansion coefficient in terms of R~ n functions [Eq. (J7J)]. Note that 
the fit of Eq. ([2]) is mathematically equivalent to evaluating the overlap integral between 
the angular functions and V(R, 9a, &b, 4>) by Gauss-Legendre quadrature. The fit was done 
using a linear least-squares procedure in which the ab initio points were weighted with 
the appropriate quadrature weights and a factor of R 3 . The /^-dependent factor is chosen 



9 



because the leading dipole-dipole interaction decays as R~ 3 . 

In principle, our long-range expansion is valid for all terms up to La — Lb — 10, with 
eight possible values of n for each (La, Lb, L) term. However, the inclusion of high powers of 
1/R may lead to unphysical results even for the low-n coefficients, which are considered the 
most important. Thus, we must carefully choose which R~ n A LmLb ^ l (Q a,$b,4>) functions to 
include in the fit. After extensive testing, we found that the best analytical fit is obtained 
for n < 14. This result is based on a thorough examination of both the stability of the 
fit, i.e. how much the Cl a ,l b ,Lm coefficients vary when adding more R~ n functions, and the 
RMS error in the data points. The final fit gave a RMS error of 4.6 xlO -8 Eh (0.03%) for 
a total of 10648 ab initio points. The RMS difference between the analytical potential and 



the size-consistency corrected long-range points of Ref. |26J], which served as test points, was 
4.8xl0~ 7 Eh (0.24%). Note that the latter error is, in part, due to the weaker convergence 



thresholds used in the calculations of Ref. 



through EPAPS [46|. 



. The Cl A: l B: l !TI fit coefficients are available 



B. Short-range S = 2 potential 

For the short range of the quintet surface, we used the size-consistency corrected 
RCCSD(T) data of Dhont et al. fl, calculated at R values from 4.0 to 16.0 ao- The 
angular grid consisted of 11 points in 9a and 8b, ranging from 0° to 180° in steps of 20° 
with an additional point at 90°. The grid in ranged from 0° to 180° in steps of 22.5°. 
The short-range potential was first expanded in terms of A LaLb;M (9a, 9b,4>) functions [Eq. 
([3])] and then transformed to Eq. ([2]). Instead of using the two-step spline-based approach 



described in Ref. 26j, we employed a weighted least squares fitting procedure to determine 
the vl A: l B: m(R) coefficients for each R. In order to perform the fit, we first calculated op- 
timal quadrature weights for the grid points in (8a, 0b), of which the details are given in 
the Appendix. We then attempted to fit the RCCSD(T) points by an expansion in terms of 
Al A) l b ,m {9 a, &b, 0) functions, weighting each point with the appropriate quadrature weights. 
High-energy points (> 0.1 Eh), which are not of practical importance in bound-state and 
scattering calculations, were excluded from the fit. It was found, however, that the least 
squares problem of Eq. is ill-conditioned for max(L J 4,L^) > 9 due to both the choice 
of grid points (the angle is undefined if 6a or 9b equals 0° or 180°) and the omission of 
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points at high energies. We therefore employed a modified fitting scheme to regularize the 
least squares problem such that all functions up to La = Lb = 10 and M = 8 could be in- 
cluded. This was done by means of a Tikhonov regularization method in which the term 
Sl a l b m \ a (^A + L B )v l a ,l b ,m (R)\ 2 was added to the residual. The factor of a(L A + L B ), 
with a = 2 x 10~ 4 , ensures that strong oscillations (associated with large La and L B ) are 
damped out in the fit. The resulting vl a ,l b ,m(R) fit coefficients were then transformed to 
v l a ,l b ,l(R) coefficients using Eq. (jSJ). Overall, this fitting procedure gave a RMS error of 



9.8xl(T 6 E h (0.21%) based on 2 
retrieved via the EPAPS system 



275 ab initio points. The vl a ,l b ,l(R) coefficients can be 



46|. 



The vl a .l b ,l{R) expansion coefficients were interpolated in R using the reproducing ker- 
nel Hilbert space (RKHS) method with the reproducing kernel for distancelike variables 



491 ] . The RKHS parameter m, which determines the power with which the interpo- 
lated function decays between the grid points, was set to the leading power in 1/R for each 
(La, Lb, L) term. For instance, the ■Um(-R) coefficient containing the electrostatic dipole- 
dipole interaction was interpolated with m — 3, the isotropic vqoo(R) term was interpolated 
with m = 6, and so on. In all cases, the RKHS smoothness parameter was set to 2. 

Finally, we matched the short-range and long-range expansions of the RCCSD(T) quintet 
potential using an independent switching function f(R) that changes smoothly from to 1 



on the interval a < R < b: 



f(R) 



if R < a 

1 if R > b (8) 
\ + \ sin ^ (3 — sin 2 ^ ) otherwise, 



with x = — b l + ( R a ) . The function is such that the first three derivatives at R = a and 

o—a 

R = b are zero. We used Eq. (jSJ) to switch the potential between a = 8 and b = 12 a . The 
total S = 2 potential energy surface may now be expressed as follows: 

V(R, 6 A , &b, <P) = [1 - f(R)]V sr (R, Q A , Ob, <P) + f(R)V lr (R, 6 A , 6 B , <P), (9) 

where V sr refers to the short-range expansion of Eq. (j2J) and V\ r to the long-range expansion 
of Eqs. © and Q. 
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C. Short-range S = 0, 1 potentials 



As already mentioned in Section [TTJ the singlet and triplet potentials were obtained 
from the quintet RCCSD(T) potential by adding the energy differences at the CASPT2 
or CASPT3 level of theory. We fitted these exchange splittings {YcASPTn ~ ^casptu) di- 
rectly in terms of Al A} l b ,l{0a, @b, 4>) functions, weighting each point with the corresponding 
Gauss-Legendre and Gauss-Chebyshev quadrature weights. In all cases, the fit error was 
largest at 4.0 ao and rapidly decreased as a function of R. For instance, the RMS errors for 
the singlet-quintet CASPT2 and CASPT3 splittings were 1.3xl0~ 3 E h (4.6%) and 1.2xl(T 3 
E h (4.7%) at 4.0 a , l.lxlO" 5 E h (0.10%) and 7.8xl0~ 6 E h (0.09%) at the neighboring grid 
point of 4.5 oq, and 2.3xl0~ 8 E h (0.007%) and 1.9xl0~ 8 E h (0.007%) near the van der 
Waals minimum at 6.5 ao- For the triplet-quintet CASPT2 and CASPT3 exchange split- 
tings, the RMS errors were 6.9xl0- 4 E h (3.2%) and 7.9xl0~ 3 E h (4.4%) at 4.0 a , 4.3xl0~ 6 
E h (0.06%) and 5.1xl0~ 6 E h (0.08%) at 4.5 a , and 2.1xl0- 8 E h (0.01%) and 6.0xl0~ 8 
Eh (0.03%) at 6.5 ao- All errors were calculated from 1331 ab initio points per R value, 
with the exception of R = 4.0 a , where three points failed to converge. The vl*,l b ,l(R) fit 
coefficients for the CASPTn energy splittings are available through EPAPS 46]. 

The vl a ,l b ,l(R) coefficients were interpolated in R using the RKHS method. For all 
(La, Lb, L) terms we set the RKHS parameter m to 14 and the smoothness parameter to 
2. The value of m = 14 ensures that all coefficients decay as R~ 15 beyond the outermost 
grid point, thus decaying faster than any of the long-range terms included in the fit of Eq. 
(|7|). In addition, we found that the interpolation with m = 14 gives the smallest RMS 



error in the ab initio points of Ref. 



261 ] . The expanded CASPTn splittings were added to 



the RCCSD(T) potential of Eq. (jUJ) to obtain the final singlet and triplet potential energy 
surfaces. 



IV. RESULTS AND DISCUSSION 



The main features of the singlet, triplet, and quintet potentials have already been de- 
scribed in Ref. 26j , and therefore we only briefly mention them here. Our S = 2 potential is 
characterized by a van der Waals minimum at R e = 6.61 a with a well depth of D e = —675 



cm 1 . It should be noted that Dhont et al. 



261 ] reported a slightly different R e value of 6.60 
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do- The minimum corresponds to a linear geometry (9a = 9b = <P = 0°) in which the two NH 



dipoles are aligned. Zuchowski et al. 



39j have recently shown that D e changes to —693 cm 1 



if the aug-cc-pVQZ basis is used and the RCCSD(T) calculations are performed without the 
frozen-core approximation. They also demonstrated from symmetry-adapted perturbation 
theory (SAPT) calculations that the main contributions to D e are the electrostatic (—899 
cm -1 ) and dispersion (—432 cm -1 ) interactions. The total SAPT exchange- repulsion energy 
at the minimum was found to be 874 cm -1 391 ] . 



The V 2 ' s=0 (V^ =0 ) and V^ =l (V^ =1 ) surfaces also exhibit a van der Waals minimum at 
9a = 9b = 4> = 0°, located at R e = 6.50 (6.51) and 6.54 (6.55) ao, respectively. These 
distances are 0.01 - 0.02 ao different from the R e values reported by Dhont et al. 26]. 
Furthermore, the singlet and triplet dimers may form the chemically stable N 2 H 2 molecule, 
which is reflected in the strongly attractive nature of these potentials at short intermolecular 
separations. The most favorable geometries for the S = and 1 states at short distances 
are found near 9 a = 9b = 90°. 



A. Long-range potential 

Before discussing the analytical fit results, we first address the size-consistency problem 
occurring at the RCCSD and RCCSD(T) levels of theory. Figure [TJ shows the isotropic part 
of the quintet potential, vooo(R), between R = 15 and 30 ao- The lack of size consistency is 
most apparent at the RCCSD level, giving rise to an error of —1.07 cm -1 at long range. The 
inclusion of triple excitations reduces the problem significantly, but in fact overcompensates 
for the RCCSD error by +0.11 cm -1 . The uncorrected isotropic RCCSD and RCCSD (T) 
potentials cross at R ~ 11 a . After subtracting the size-consistency errors from all ab 
initio points, both the RCCSD and RCCSD (T) potentials smoothly converge to zero at 
long range. It can also be seen that these corrected data are in very good agreement with 
the corresponding spin- unrestricted CC results at R — 15 and 30 a . 

The main fit results for the (size-consistency corrected) RCCSD (T) long-range potential 
are presented in Table [TO A total number of 588 C La,Lb ,L,n coefficients was included in 
the long-range fit (La, Lb < 10 and n < 14), but here we list only the most important 
terms. Table [III also shows the results obtained from first and second-order perturbation 
theory (PT). It can be seen that the fitted electrostatic terms agree very well with the PT 
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coefficients, in particular with the data calculated at the PT-RCCSD(T) level of theory. 
For the induction and dispersion terms we find some significant discrepancies, but the most 
important second-order fit coefficients (those with n = 6) show satisfactory agreement with 
PT-DFT. It should be noted that, for the fitted coefficients, no distinction can be made 
between induction and dispersion contributions. For the isotropic C 0i o,o,6 term, the PT- 
DFT calculations give a dispersion coefficient of 39.86 a.u. and an induction term of 6.99 
a.u. 

As an indication of the relative importance of the different C£ Ai L B ,L,n coefficients, we 
explicitly give their contributions to the potential at R = 30 a (see Table HI])- These 
contributions, VL A ,L B ,L,n(R), were calculated as follows: 

I C^ u I 

v La , Lb ,UR) = n La>l J LA '^' L ' n \ (io) 

where N La>Lb;L = [^/{2L A + 1){2L B + 1)(2L + l)] 1 / 2 is the norm of the angular 
Al a ,l b ,l(0a, 8b, 0) functions. It is clear that the n = 3 dipole-dipole interaction domi- 
nates the potential by at least one order of magnitude, followed by the electrostatic dipole- 
quadrupole term. The main second order term is the isotropic n = 6 interaction, which, 
at 30 do, is still larger than the electrostatic n = 5 terms. The fact that the fitted 61^2,3 
and Co,o,o,6 coefficients give the largest contributions in first and second order, respectively, 
indicates that the fit is not only numerical, but also physically meaningful. Thus, we may 
safely extrapolate the potential from 30 ao to larger R values. 

Figure [2] shows the independence of the fitted RCCSD(T) potential for two specific 
orientations {8 a, 8 b, </>)• For the linear geometry, with 8 a — 8 B — <fi — 0°, the leading dipole- 
dipole interaction is maximally attractive, while for 8a = 9b — <P = 90° the dipole-dipole 
term is zero. It can be seen that the 6*1^2,3 coefficient dominates the long-range potential 
beyond R ~ 12 a . Figure [2] also compares the total long-range expansion with the ab initio 
data, illustrating the region of validity of Eq. ([7]). It should be noted that, on the scale 
of the figure, the short-range expansion of Eq. (T5]) is indistinguishable from the total fitted 
potential of Eq. ([9]), and thus the short-range expansion is not explicitly shown. The long- 
range fit is very accurate for intermolecular distances larger than 8 a , which suggests that 
short-range (exchange and charge penetration) effects are only significant for R < 8 ao- This 
also justifies our choice of switching the potential from the short-range to the long-range 
expansion between 8 and 12 a . 
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B. Short-range potentials 



Although the 5 = 0, 
reported by Dhont et al. 



, and 2 potentials obtained in this work are very similar to those 



26| , there are some notable differences at very short intermolecular 



distances. The differences are most pronounced at R = 4.0 a , where the potentials exhibit 
the highest anisotropy. Figure [3] compares the two fit results for the quintet state as a 
function of 9a and 9b, with R = 4.0 a and <fi = 0°. Note that both surfaces were obtained 
from the same set of ab initio data. The fit of Ref. [26| shows more oscillatory behavior than 
our present result, in particular near (9a, Ob) = (180°, 150°) and (30°, 0°). Furthermore, the 
potential of Dhont et al. has a local maximum around (150°, 30°) that is clearly unphysical 
in nature. Similar patterns are found for the triplet and singlet states, as can be seen in Figs. 
H] and [51 The 5 = and 1 potentials of Ref. [26| exhibit more pronounced oscillations and 
local maxima, indicating more unphysical behavior. We therefore conclude that, in addition 
to the more accurate long-range potential, the fit of the short-range NH-NH potentials is 
also improved in the present work. 



V. CONCLUSIONS AND OUTLOOK 



We have constructed four-dimensional potential energy surfaces for the singlet, triplet, 
and quintet states of NH( 3 £~) - NH( 3 S~) based on high-level ab initio calculations. All 
potentials were fitted in terms of coupled spherical harmonics in the angular coordinates, and 
the long range was further expanded as a power series in 1/R. Prior to fitting, the ab initio 
data were corrected for a size-consistency error of 0.5 x 10~ 6 Eh occurring at the RCCSD(T) 
level of theory. The fitted long-range coefficients were found to be in good agreement with 
the results obtained from first and second-order perturbation theory. 

Future work is planned to study the evaporative cooling process of NH, which requires 
knowledge of elastic and inelastic cross sections at (ultra) low temperatures. Rate constants 
and cross sections for (cold) reactive NH + NH collisions will also be calculated. Finally, we 
aim to explore the possibilities of cold controlled chemistry by investigating the influence of 
external fields. 
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APPENDIX 

In this Appendix, we describe how we optimized the quadrature weights Wi for the 
integration of Legendre polynomials Pi(x) on a given grid of mutually distinct points Xi 
(i = l,...,n): 

/l n 
P l {x)dx = 25 l>0 nJ2w i P l {x i ). (A.l) 
1 i=i 

We define the optimization as a minimization of the sum of square residuals \r\: 

\r\ = \Aw - b\, (A.2) 

where A is an (l max + 1) x n matrix with elements An = P\{xi) (/ = 0, . . . , l max ), w is a 
vector of length n containing the quadrature weights Wi, and 6 is a vector of length l max + 1 
with elements bi = 2^ - m the case of an n-point Gauss-Legendre quadrature, x» and W{ 
are chosen in such a way that the integration is exact, i.e., \r\ = 0, for all polynomials up to 
degree l max = 2n — l. For arbitrary, mutually distinct points x«, we may calculate the weights 



as w = A _1 6, since A is regular for l max = n — 1 (see p. 145 of Ref. [32]). This results in a 
quadrature that is exact up to (at least) degree n — 1. Instead of using a quadrature that is 
exact for l max = n — 1 and most likely unsuitable for higher degree polynomials, we choose a 
compromise quadrature that is reasonable for l max > n — 1 at the expense of not being exact 
for lower degree polynomials. This may be achieved by linear least squares minimization 
of \r\. However, we prefer to use a quadrature that is exact for constant functions (1 = 0), 
which requires a minimization of \r\ with the constraint that Y17=i w i = 2- For this purpose 
we take 

w = Wq + (A. 3) 
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with (wo)i = 2/n for all % = l,...,n and Y^i=i( Wl )i = ®- This may be rewritten as 
WqW± = 0, with Wq denoting the transpose of Wq. We can now expand w± in an orthogonal 
basis {qi, i = 2, . . . , n} of vectors q^ that are perpendicular to w : 

n 

w ± = ^2 q^ = Qc. (A.4) 

i=2 

We observe that the first row of the matrix A is proportional to Wq, and thus the vectors 
qi can be generated by Gram-Schmidt QR-factorization of A T : 

A T = QR. (A.5) 

Here, Q is an n x n ortho normal matrix and R is an n x (l max + 1) upper triangular 
matrix. The columns i — 2, . . . , n of Q form the matrix Q of Eq. (IA.4I) . In order to find 
the expansion coefficients c, we now remove the first row of A and the first element of b, 
yielding the (l max x n) matrix A and the null vector b of length l ma x, respectively, and define 
the residual r = Aw. Substitution of Eq. flA.31) gives 

\r\ = \Aw + AQc\, (A. 6) 

which can be minimized in a standard least squares procedure to obtain the expansion 
coefficients c. Finally, substitution of Eq. flA.41) into flA.31) gives the total optimal quadrature 
weights. In the present work, we have employed this method to generate optimal weights 
for the short-range quintet potential with n — 11 and l max = 16. 
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TABLE I: Size-consistency errors (AE) for the NH-NH system at the RCCSD and RCCSD(T) 
levels of theory. The basis sets correspond to the (aug)-cc-pVnZ (n = double, triple, quadruple, 



quintuple) sets of Dunning [3l|. The errors are calculated as the difference between the energy of 
the separate monomers and the energy of the supersystem NH- • • NH at an intermolecular distance 
of 30000 a . All values are in 10~ 6 E h . 



Basis set 



AE RCCSD 



AE RCCSD(T) 



cc-pVDZ 
cc-pVTZ 
cc-pVQZ 
cc-pV5Z 



-3.15067 
-4.25041 
-4.70853 
-4.92130 



-0.50946 
-0.01069 
0.36976 
0.62672 



aug-cc-pVDZ 
aug-cc-pVTZ 
aug-cc-pVQZ 
aug-cc-pV5Z 



-4.04159 
-4.88230 
-5.01375 
-5.03493 



0.01944 
0.51290 
0.68981 
0.75827 
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TABLE II: Most important long-range coefficients obtained from the fit and from perturbation 
theory, and their contributions at 30 cto- The order of importance is based on the value of n, and 
for each n only the four largest terms are given. Terms labeled with an asterisk are first-order 
(electrostatic) interactions. All values are in atomic units. Numbers in parentheses denote powers 



of 10. 





L B 


L 


n 




r PT~RCCSD(T) 


sjPT—DFT 
L A .Lfi.L.yi 


V LA ,L B ,L,n( m °o) 


1 


1 


2 


3* 


1.9697(+0) 


1.9715(+0) 


2.0127(+0) 


3.8551(-05) 


1 


2 


3 


4* 


-2.8394(+0) 


-2.8597(+0) 


-3.0642(+0) 


1.2127(-06) 


1 


3 


4 


5* 


1.6637(+1) 


1.6761(+1) 


1.7103(+1) 


1.7654(-07) 


2 


2 


4 


5* 


-5.6953(+0) 


-5.4312(+0) 


-6.1080(+0) 


5.5389(-08) 











6 


4.7270(+l) 




4.6852(+l) 


2.2986(-07) 


1 


4 


5 


6* 


-5.4131(+1) 


-5.5049(+l) 


-5.7422(+l) 


1.5274(-08) 





2 


2 


6 


1.2657(+1) 




1.5681(+1) 


1.2309(-08) 


2 


3 


5 


6* 


3.6904(+l) 


3.9347(+l) 


4.2140(+1) 


9.1458(-09) 





1 


1 


7 


-1.8433(+2) 




-8.2153(+1) 


9.9596(-09) 


1 


2 


3 


7 


-3.4979(+2) 




1.5651(+1) 


5.5331(-09) 





3 


3 


7 


-1.0784(+2) 




-7.9522(+l) 


2.4971(-09) 


3 


3 


6 


7* 


3.1701(+2) 


3.3946(+2) 


3.4622(+2) 


2.0359(-09) 











8 


9.2546(+2) 




1.1077(+3) 


5.0003(-09) 





2 


2 


8 


3.9371(+3) 




1.4208(+3) 


4.2544(-09) 


1 


1 


2 


8 


4.5792(+3) 




-1.0618(+2) 


3.6882(-09) 


2 


2 


4 


8 


-4.4826(+3) 




6.2384(+2) 


1.6146(-09) 





1 


1 


9 


3.0500(+4) 




-3.1644(+3) 


1.8310(-09) 


1 


2 


3 


9 


9.9936(+4) 




2.2929(+3) 


1.7565(-09) 


1 


2 


1 


9 


-2.3093(+4) 




-6.1280(+2) 


6.2000(-10) 





3 


3 


9 


9.0400(+3) 




-5.6295(+3) 


2.3259(-10) 
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FIG. 1: Isotropic part of the quintet potential calculated at the CCSD and CCSD(T) levels of the- 
ory. The data labeled "RCCSD" and "RCCSD(T)" correspond to the uncorrected spin-restricted 
data, "RCCSD-AE" and "RCCSD(T)-AE" to the size-consistency corrected data, and "UCCSD" 
and "UCCSD(T)" to the spin-unrestricted results. 

FIG. 2: independent quintet potential for two selected orientations (0j^,8b,(P). The solid lines 
correspond to the total fitted potential, the dashed lines to the long-range potential, and the 
dotted lines to the long-range dipole-dipole interaction. 

FIG. 3: (Color) Cuts of the quintet potential (in cm -1 ) for R = 4.0 ag and cj) = 0°. The left panel 



shows the fit obtained in this work and the right panel shows the results of Ref. 



26]. 



FIG. 4: (Color) Cuts of the triplet potential (in cm -1 ) for R = 4.0 ao and (j) = 0°, calculated using 
Eq. ([1]). Th e up per panels correspond to the present work and the lower panels to the work of 
Dhont et al. 

FIG. 5: (Color) Cuts of the singlet potential (in cm -1 ) for R = 4.0 ao and <p = 0°, calculated using 
Eq. ([1]). Th e up per panels correspond to the present work and the lower panels to the work of 
Dhont et al. 12611 . 
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Janssen et al. Fig. [TJ 
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Janssen et al. Fig. [2J 
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Janssen et al. Fig. [3J 
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Janssen et al. Fig. HI 
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Janssen et al. Fig. [5j 



28 



